#ifndef __RANDVALS_DOT_HPP_GUARD__
#define __RANDVALS_DOT_HPP_GUARD__
////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////
#include <stdint.h>
////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////
#include <boost/random/mersenne_twister.hpp>
#include <boost/random/uniform_int.hpp>
#include <boost/random/uniform_real.hpp>
#include <boost/random/variate_generator.hpp>
#include <boost/math/distributions/pareto.hpp>
//#include <boost/random.hpp>
////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////
uint64_t RandU64( void );
////////////////////////////////////////////////////////////////////////////////
////////////////////////////////////////////////////////////////////////////////
class RandomNumberGenerator {
    
    const uint64_t seed;
    boost::mt19937                  generator;
    boost::uniform_int< uint64_t >  uniIntDist;
    boost::uniform_real< double  >  uniDoubleDist;
    boost::variate_generator< boost::mt19937, boost::uniform_int< uint64_t > > randomU64;
    boost::variate_generator< boost::mt19937, boost::uniform_real< double > >  randomDouble;
    
public:
    RandomNumberGenerator( void ) :
    seed( time( NULL ) ),
    generator( seed ),
    uniIntDist( 0ULL, 0xFFFFFFFFFFFFFFFFULL ),
    uniDoubleDist( 0.0, 1.0 ),
    randomU64( generator, uniIntDist ),
    randomDouble( generator, uniDoubleDist )
    {}

    inline uint64_t RandU64( void ) {
        return randomU64();
    }
    inline double   RandDouble( void ) {
        return randomDouble();
    }
};

class PowerLaw {
    boost::math::pareto_distribution<> dist;
    RandomNumberGenerator rng;
public:
    PowerLaw( double a ) : dist( 1.0, a ) {
        assert( a > 0.0 );
    }
    double Sample( void ) {
        const double p = rng.RandDouble();
        const double x = boost::math::quantile( dist, p );
        return x;
    }
    uint64_t IntSample( void ) {
        double x = Sample();
        x = floor( x + 0.5 );
        assert( x > 0ULL );
        return uint64_t( x );
    }
};


#endif